
#include <mecafluide.hpp>
#include <feel/feelfilters/geotool.hpp>
#include <iostream>

//#include <feel/feelfilters/gmshhypercubedomain.hpp>
#include <feel/feelfilters/gmshsimplexdomain.hpp>

#include <feel/feeldiscr/operatorinterpolation.hpp>

#include <feel/feeldiscr/operatorlagrangep1.hpp>

//#include <mecafluide_cl.cpp>

#include <research/life-workspace/Fluid-Structure/functionSup.cpp>

#include <research/life-workspace/Fluid-Structure/methodsnum.hpp>


namespace Feel {

namespace FSI {

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    NavierStokes<Dim,Order,Entity>::NavierStokes( application_ptrtype __appli )
        :
        M_lifeApplication(__appli),
        M_backend( backend_type::build( this->application()->vm() ) ),
        meshSize( this->application()->vm()["fluide-hsize"].template as<double>() ),
        mesh_fluid( this->createFluidMesh( meshSize,"domainFluid" )),
        mesh_reference_fluid( this->createFluidMesh( meshSize,"domainFluidRef" )),
        //-----------------------------------------------------------------------------//
        //M_mesh_interface(this->createInterfaceMesh(meshSize, "InterfaceFluid" ) ),
        M_Xh_interface(new space_fluid_velocity_type( mesh_fluid )),
        M_velocity_interface(new element_fluid_velocity_type( M_Xh_interface, "velocity_interface" ) ),
        M_Xh_interface_ref(new space_fluid_velocity_type( mesh_reference_fluid )),
        M_velocity_interface_ref(new element_fluid_velocity_type( M_Xh_interface_ref, "velocity_interface_ref" ) ),
        //-----------------------------------------------------------------------------//
        mesh_struct_top(this->createStructMesh( meshSize,0.5/*1.0*/, "boundaryMoveH" ) ),
        mesh_struct_bottom(this->createStructMesh( meshSize,-0.5/*0.0*/, "boundaryMoveB" ) ),
        Sh_top(new structure_space_type( mesh_struct_top )),
        Sh_bottom(new structure_space_type( mesh_struct_bottom )),
        // reference boundary description
        M_bc_reference_top(new element_structure_type( Sh_top, "bc_reference_top" ) ),
        M_bc_reference_bottom(new element_structure_type( Sh_bottom, "bc_reference_bottom" ) ),

        M_Xh(new space_fluid_type( mesh_fluid ) ),
        M_Xh_ref(new space_fluid_type( mesh_reference_fluid ) ),
        U_fluid( new element_fluid_type(M_Xh,"U")),
        M_Ufluid_ref( new element_fluid_type(M_Xh_ref,"Uref")),
        M_identity_ale(new element_fluid_velocity_type( M_Xh->template functionSpace<0>(), "identity_ale" ) ),
        M_mesh_velocity(new element_fluid_velocity_type( M_Xh->template functionSpace<0>(), "mesh_velocity" ) ),
        M_bdf_fluid ( new bdf_type( this->application()->vm(),M_Xh,"fluid" )),
        M_depl_top(new element_structure_type( Sh_top, "bc_top" )),
        M_depl_bottom(new element_structure_type( Sh_bottom, "bc_bottom" )),
        M_aleFactory( new ALE< convex_type >( std::make_pair(0,6), mesh_reference_fluid, this->application()->vm() ) ),
        M_aux(new ale_map_element_type( M_aleFactory->getDisplacement().functionSpace(), "new ale map" ) ),
        M_bdf_ale_map ( new bdf_ale_map_type(this->application()->vm(),
                                             M_aleFactory->getDisplacement().functionSpace(),
                                             "ale_map") ),
        M_bdf_ale_identity ( new bdf_ale_velocity_type( this->application()->vm(),
                                                        M_Xh->template functionSpace<0>(),
                                                        "velocity fluid" ) ),
        M_exporter( Exporter<mesh_type>::New( /*"gmsh"*/this->application()->vm(), /*this->application()->about().appName()*/"ExportFluid" ) ),
        M_exporter_ref( Exporter<mesh_type>::New( /*"gmsh"*/this->application()->vm(), /*this->application()->about().appName()*/"ExportFluid_ref" ) )
        //M_exporter_interface(Exporter</*mesh_structure_type*/Mesh< Simplex<2,1,2> > >::New( this->application()->vm(), this->application()->about().appName()+"_interface" ) )
    {
        M_time = __appli->vm()["bdf-time-initial"].template as<double>();
        M_cpt_export=0;
        //Masse volumique
        M_rho = this->application()->vm()["rho"].template as<double>();
        M_mu = this->application()->vm()["mu"].template as<value_type>();
        M_nu = M_mu/M_rho;
        penalbc = this->application()->vm()["bccoeff"].template as<value_type>();
        Vitesse_Entree = this->application()->vm()["entree-v"].template as<value_type>();
        Pression_Sortie = this->application()->vm()["sortie-p"].template as<value_type>();
        M_weakCL = this->application()->vm()["weakCL"].template as<bool>();
        M_isStat = this->application()->vm()["stationnaire"].template as<bool>();
        M_isMoveDomain = this->application()->vm()["movedomain"].template as<bool>();
        M_isFixeBoundary = false;

        *M_bc_reference_top = vf::project( Sh_top, elements(Sh_top->mesh()), vf::constant(1.0) );
        *M_bc_reference_bottom = vf::project( Sh_bottom, elements(Sh_bottom->mesh()), vf::constant(0.0) );
        M_referencePolyBoundarySet.push_back(*M_bc_reference_bottom);
        M_referencePolyBoundarySet.push_back(*M_bc_reference_top);


        M_LookComment = this->application()->vm()["fluide-verbose"].template as<bool>();//false;
#if 0
        // set up the non linear solver
        M_backend->nlSolver()->residual = boost::bind( &NavierStokes<Dim, Order, Entity>::updateResidual,
                                                       boost::ref( *this ), _1, _2 );
        M_backend->nlSolver()->jacobian = boost::bind( &NavierStokes<Dim, Order, Entity>::updateJacobian,
                                                       boost::ref( *this ), _1, _2 );

#endif
        std::cout << "constructeur NS finsih\n";

    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    typename NavierStokes<Dim,Order,Entity>::mesh_ptrtype
    NavierStokes<Dim,Order,Entity>::createFluidMesh( double meshSize,std::string name )
    {
        if (M_LookComment) std::cout << "[NavierStokes] : createFluidMesh start\n";

        this->changeRepository();


        //Geometry
#if 0
        GeoTool::Node x1(0,0);
        GeoTool::Node x2(6,1);
#else
        GeoTool::Node x1(0,-0.5);
        GeoTool::Node x2(6,0.5);

#endif
        //GeoTool::Node x1(-0.5,-0.5);
        //GeoTool::Node x2(1,1.5);


        GeoTool::Rectangle R( meshSize,"LEFT",x1,x2);
        R.setMarker(_type="line",_name="Entree",_marker4=true);
        //R.setMarker(_type="line",_name="Paroi",_marker1=true,_marker3=true);
        R.setMarker(_type="line",_name=/*"ParoiB"*/"Paroi",_marker1=true);
        R.setMarker(_type="line",_name=/*"ParoiH"*/"Paroi",_marker3=true);
        R.setMarker(_type="line",_name="Sortie",_marker2=true);
        R.setMarker(_type="surface",_name="OmegaFluide",_markerAll=true);

        GeoTool::Node x_centre(2.5,0.5);
        GeoTool::Node x_bord(2.65,0.5);
        GeoTool::Circle C( meshSize/3.,"LEFT2",x_centre,x_bord);
        C.setMarker(_type="line",_name="Paroi",_markerAll=true);
        C.setMarker(_type="surface",_name="OmegaFluide",_markerAll=true);

        GeoTool::Node x_centre2(2.5,0.3);
        GeoTool::Node x_bord2(2.7,0.3);
        GeoTool::Circle C2( meshSize/3.,"disque2",x_centre2,x_bord2);
        C2.setMarker(_type="line",_name="Paroi",_markerAll=true);
        C2.setMarker(_type="surface",_name="OmegaFluide",_markerAll=true);

        if (M_LookComment) std::cout << "[NavierStokes] : createFluidMesh finish\n";

        //return (R-C).createMesh<mesh_type>(name);
        return R.createMesh<mesh_type>(name);

    } // NavierStokesStokes::createFluidMesh


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    typename NavierStokes<Dim,Order,Entity>::mesh_structure_ptrtype
    NavierStokes<Dim,Order,Entity>::createInterfaceMesh( double meshSize, std::string name )
    {

        if (M_LookComment) std::cout << "[NavierStokes] : createInterfaceMesh start\n";

        this->changeRepository();

        GeoTool::Node x1(0,1);
        GeoTool::Node x2(6,1);
        GeoTool::Line G1( meshSize,"interface_haut",x1,x2);

        GeoTool::Node x3(0,0);
        GeoTool::Node x4(6,0);
        GeoTool::Line G2( meshSize,"interface_bas",x3,x4);

        if (M_LookComment) std::cout << "[NavierStokes] : createInterfaceMesh finish\n";

        return (G1+G2).createMesh<mesh_structure_type>(name);

    }

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    typename NavierStokes<Dim,Order,Entity>::mesh_structure_ptrtype
    NavierStokes<Dim,Order,Entity>::createStructMesh( double meshSize,double __ypos, std::string name )
    {
        if (M_LookComment) std::cout << "[NavierStokes] : createStructMesh start\n";

        this->changeRepository();

        GeoTool::Node x1(0,__ypos);
        GeoTool::Node x2(6,__ypos);
        GeoTool::Line R( meshSize,"kokok",x1,x2);

        if (M_LookComment) std::cout << "[NavierStokes] : createStructMesh finish\n";

        return R.createMesh<mesh_structure_type>(name);
    }



    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::exportResults()
    {
        if (M_LookComment) std::cout << "[NavierStokes] : Export NS start\n";

        this->changeRepository();
        //{ this->changeRepository("olaExport");M_cpt_export=0;}


#if 0
#if 0
        MeshHighOrder< /*convex_type*/Simplex<Dim,1,Dim> > ho_mesh ( mesh_fluid );
        auto ho_fluid_mesh = ho_mesh.getMesh();

        typedef bases<Lagrange<Order, Vectorial, Continuous, PointSetFekete> > pN_ale_basis_type;
        typedef FunctionSpace< mesh_type, pN_ale_basis_type, double> pN_visualize_functionspace_type;
        typedef boost::shared_ptr<pN_visualize_functionspace_type> pN_visualize_functionspace_ptrtype;
        typedef typename pN_visualize_functionspace_type::element_type pN_element_type;

        pN_visualize_functionspace_ptrtype visH = pN_visualize_functionspace_type::New( ho_fluid_mesh );
        pN_element_type aux_element( visH, "aux");
        aux_element = vf::project( visH, elements(visH->mesh()), vf::idv(U_fluid->template element<0>()) );
#endif
        //typedef space_fluid_ptrtype functionspace_ptrtype;
        typedef space_fluid_velocity_ptrtype functionspace_ptrtype;
        //typedef space_fluid_type functionspace_type;
        typedef space_fluid_velocity_type functionspace_type;
        functionspace_ptrtype Wh = functionspace_type::New( /*aux_element.functionSpace()->mesh()*/mesh_fluid );
        //functionspace_ptrtype Wh = functionspace_type::New( aux_element.functionSpace()->mesh() );
        //space_fluid_ptrtype Wh = space_fluid_type::New( aux_element.functionSpace()->mesh() );
        OperatorLagrangeP1<functionspace_type> I( Wh, M_backend );

        typename OperatorLagrangeP1<functionspace_type>::dual_image_space_ptrtype Yh( I.dualImageSpace() );
        //exporter->step(tn)->setMesh( Yh->mesh() );

        M_exporter->step( M_cpt_export )->setMesh(Yh->mesh() );
        M_exporter->step( M_cpt_export )->add( "vitesse", U_fluid->template element<0>() );
        M_exporter->save();

        space_fluid_pressure_ptrtype Wh2 = space_fluid_pressure_type::New( mesh_fluid );
        OperatorLagrangeP1<space_fluid_pressure_type> I2( Wh2, M_backend );
        typename OperatorLagrangeP1<space_fluid_pressure_type>::dual_image_space_ptrtype Yh2( I2.dualImageSpace() );
        M_exporter_ref->step( M_cpt_export )->setMesh( Yh2->mesh() );
        M_exporter_ref->step( M_cpt_export )->add( "pression", U_fluid->template element<1>() );
        M_exporter_ref->save();

#else
        M_exporter->step( M_cpt_export )->setMesh( mesh_fluid );
        M_exporter->step( M_cpt_export )->add( "vitesse", U_fluid->template element<0>() );
        M_exporter->step( M_cpt_export )->add( "pression", U_fluid->template element<1>() );
        M_exporter->step( M_cpt_export )->add( "mesh_velocity", *M_mesh_velocity );
        M_exporter->step( M_cpt_export )->add( "aux", *M_aux );
        M_exporter->step( M_cpt_export )->add( "mesh_velocity_interf", *M_velocity_interface );
        M_exporter->save();
#endif




#if 1
        //M_exporter_ref->step( M_cpt_export )->setMesh( mesh_reference_fluid );
        M_exporter_ref->step( M_cpt_export )->setMesh( M_Ufluid_ref->mesh() );
        M_exporter_ref->step( M_cpt_export )->add( "u", M_Ufluid_ref->template element<0>() );
        M_exporter_ref->save();
#endif

#if 0
        M_exporter_interface->step( M_cpt_export )->setMesh( M_mesh_interface);
        M_exporter_interface->step( M_cpt_export )->add( "u", *M_velocity_interface);
        M_exporter_interface->save();
#endif

        ++M_cpt_export;

        if (M_LookComment) std::cout << "[NavierStokes] : Export NS finish\n";

    } // NavierStokes::export


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::changeRepository(std::string pathLoc)
    {
#if 0
        this->application()->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/"+ pathLoc )
                                               % this->application()->about().appName()
                                               % convex_type::name()
                                               % Order
                                               % this->meshSize
                                               );

#endif
    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    double
    NavierStokes<Dim, Order, Entity>::diffL2( element_fluid_type& U1, element_fluid_type& U2 )
    {
        using namespace vf;

        element_fluid_0_type u1 = U1.template element<0>();
        element_fluid_1_type p1 = U1.template element<1>();
        element_fluid_0_type u2 = U2.template element<0>();
        element_fluid_1_type p2 = U2.template element<1>();

        mesh_ptrtype mesh = U1.functionSpace()->mesh();

        double __norm = math::sqrt( integrate( elements(mesh),// _Q<2*(Order+OrderGeo-1)>(),
                                               trans(idv(u1))*idv(u1) +idv(p1)*idv(p1) ).evaluate()(0,0) );
        double __diff = math::sqrt( integrate(elements(mesh),// _Q<2*(Order+OrderGeo-1)>(),
                                              trans(idv(u1)-idv(u2))*(idv(u1)-idv(u2))
                                              +(idv(p1)-idv(p2))*(idv(p1)-idv(p2)) ).evaluate()(0,0) );
        if ( __norm != 0 ) __diff /= __norm;

        return __diff;
    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::solveFluid(std::string __PDEtype)
    {

        if (M_LookComment) std::cout <<"[NavierStokes] : exec start \n";

        this->changeRepository();

        // vector solution of last step time
        vector_ptrtype Uvec( M_backend->newVector( M_Xh ) );
        *Uvec = *U_fluid;

        auto clDef = MECAFLUIDE_CL(M_time);

        MethodsNum< self_ptrtype > meth(this->shared_from_this());

        //std::string PDEtype = "Stokes";
        if (__PDEtype == "Stokes")
            meth.linearSolver(clDef,Uvec);
        else if (__PDEtype == "Navier-Stokes")
            meth.AlgoNewton(clDef,Uvec);

        //meth.AlgoNewton2(clDef,Uvec);

        *U_fluid = *Uvec;

        if (M_LookComment) std::cout <<"[NavierStokes] : exec finsh \n";

    }

    //---------------------------------------------------------------------------------------------------------//

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateLinearPDE(const vector_ptrtype& X,sparse_matrix_ptrtype& A , vector_ptrtype& F)
    {
        updateStokes(X,A,F);
    }

    //---------------------------------------------------------------------------------------------------------//
    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateALEmap()
    {

        using namespace Feel::vf;

        if (M_LookComment) std::cout <<"[NavierStokes] : update ALE map start \n";

        this->changeRepository();


        //structure_space_ptrtype Sh_top = M_depl_top->functionSpace();
        //structure_space_ptrtype Sh_bottom = M_depl_bottom->functionSpace();

        //Attention diferent celon que l'on donne le deplacement ou
        //la position de l'interface
        element_structure_type bc_disp_top( Sh_top, "bc_top" );
        //bc_disp_top = *M_depl_top;
        //bc_disp_top -= *M_bc_reference_top;
        bc_disp_top = *M_depl_top;
        bc_disp_top.updateGlobalValues();
        element_structure_type bc_disp_bottom( Sh_bottom, "bc_bottom" );
        bc_disp_bottom = *M_depl_bottom;
        //bc_disp_bottom -= *M_bc_reference_bottom;//vaut 0
        bc_disp_bottom.updateGlobalValues();
        std::vector<element_structure_type> polyBoundarySet;
        polyBoundarySet.push_back(bc_disp_bottom);
        polyBoundarySet.push_back(bc_disp_top);

        //define set of flags
        std::map< std::string, std::vector<flag_type> > flagSet;
        flagSet["fixed_bc"].push_back(mesh_fluid->markerName("Sortie"));//1
        flagSet["fixed_bc"].push_back(mesh_fluid->markerName("Entree"));//3
        flagSet["moving_bc"].push_back(mesh_fluid->markerName("ParoiB"));//2
        flagSet["moving_bc"].push_back(mesh_fluid->markerName("ParoiH"));//4

        // generate alemap with boundary on reference and actual mesh
        M_aleFactory->generateHighOrderMap( flagSet["moving_bc"], M_referencePolyBoundarySet, polyBoundarySet );

        *M_aux = M_aleFactory->getDisplacement();
        ale_map_element_type displMesh( M_aux->functionSpace() );
        displMesh = *M_aux;
        //essaye de comprendre pourquoi? !!!!!!!!!!!!!!!!!!!
        bool remesh = false;
        if ( remesh )
            {
                mesh_fluid.reset(new mesh_type());
                mesh_fluid = this->createFluidMesh( meshSize,"domainFluid" );
            }
        else
            {
                //*M_aux -= M_bdf_ale_map->unknown(0);
                displMesh -= M_bdf_ale_map->unknown(0);
            }

        //M_aux->updateGlobalValues();

        //deplacement du maillage
        M_mesh_mover.apply(mesh_fluid, /**M_aux*/displMesh );

        if ( remesh )
            {
                M_Xh = space_fluid_type::New( mesh_fluid );
                U_fluid.reset(new element_fluid_type(M_Xh, "ub"));
            }

        M_identity_ale.reset(new element_fluid_velocity_type(M_Xh->template functionSpace<0>(), "identity_ale" ) );
        *M_identity_ale = vf::project( M_Xh->template functionSpace<0>(), elements( M_Xh->mesh() ), P() );
        //M_identity_ale->updateGlobalValues();

        //std::cout << "\n vitesse mesh = " << M_bdf_ale_identity->polyDerivCoefficient(0) << "\n";

        *M_mesh_velocity = *M_identity_ale;
        //mesh_velocity.scale( M_bdf_ale_identity->polyDerivCoefficient(0) );
        //mesh_velocity.add( -1.0, M_bdf_ale_identity->polyDeriv() );
        *M_mesh_velocity -= M_bdf_ale_identity->unknown(0);
        M_mesh_velocity->scale( M_bdf_ale_identity->polyDerivCoefficient(0) );
        M_mesh_velocity->updateGlobalValues();

        //---------------------------------------------------------------

        //M_velocity_interface->container() = M_mesh_velocity->container();
        *M_velocity_interface = *M_mesh_velocity;

        //---------------------------------------------------------------

        if (M_LookComment) std::cout <<"[NavierStokes] : update ALE map finish \n";

    }

    //---------------------------------------------------------------------------------------------------------//

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::init()
    {
        if (M_LookComment) std::cout<<"[NavierStokes] : init start\n";

        this->changeRepository();

        //deplacement initial
        *M_depl_top = vf::project( Sh_top, elements(Sh_top->mesh()), vf::cst(0.0) );
        *M_depl_bottom = vf::project( Sh_bottom, elements(Sh_bottom->mesh()), vf::cst(0.0) );

        if (M_cpt_export>0) {
            // m.a.j de la map ALE
            this->updateALEmap();
        }

        // initialisation des bdfs
        *M_identity_ale = vf::project( M_Xh->template functionSpace<0>(), elements( M_Xh->mesh() ), vf::P() );
        M_identity_ale->updateGlobalValues();

        M_bdf_ale_identity->shiftRight( *M_identity_ale );
        M_bdf_ale_map->shiftRight( M_aleFactory->getDisplacement() );
        M_bdf_fluid->shiftRight( *U_fluid );

        M_bdf_ale_map->initialize( M_aleFactory->getDisplacement() );
        M_bdf_ale_identity->initialize( *M_identity_ale );
        M_bdf_fluid->initialize( *U_fluid );

        M_bdf_ale_map->start();
        M_bdf_ale_identity->start();
        M_bdf_fluid->start();

        if (M_LookComment) std::cout<<"[NavierStokes] : init finish\n";

    }


    //---------------------------------------------------------------------------------------------------------//



    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateDeplacement(element_structure_ptrtype __depl_top,
                                                        element_structure_ptrtype __depl_bottom )
    {
        if (M_LookComment) std::cout<<"[NavierStokes] : updateDeplacement start\n";

        this->changeRepository();

        M_depl_top = __depl_top;
        M_depl_bottom = __depl_bottom;
        M_depl_top->updateGlobalValues();
        M_depl_bottom->updateGlobalValues();

        if (M_isMoveDomain)
            {
                this->updateALEmap();
            }

        if (M_LookComment) std::cout<<"[NavierStokes] : updateDeplacement finish\n";
    }



    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateTime(double __t)
    {
        M_time = __t;
    }

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateBdf()
    {
        if (M_LookComment) std::cout <<"[NavierStokes] : updateBDF start\n";

        this->changeRepository();

        M_bdf_ale_identity->shiftRight( *M_identity_ale );
        M_bdf_ale_map->shiftRight( /*M_aleFactory->getDisplacement()*/*M_aux );
        //M_bdf_ale_map->shiftRight( M_aleFactory->getDisplacement() );
        M_bdf_fluid->shiftRight( *U_fluid );

        M_bdf_ale_identity->next();
        M_bdf_ale_map->next();
        M_bdf_fluid->next();

        if (M_LookComment) std::cout << "[NavierStokes] : updateBDF finish\n";
    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    typename NavierStokes<Dim, Order, Entity>::element_fluid_ptrtype
    NavierStokes<Dim, Order, Entity>::getNormalStress()
    {

        using namespace Feel::vf;

        if (M_LookComment) std::cout << "[NavierStokes] : getNormalStress start\n";
#if 0
        auto solFluid = getSolution();

        element_fluid_0_type u = solFluid->template element<0>();
        element_fluid_1_type p = solFluid->template element<1>();

        //Tenseur des deformations
        auto defv = 0.5*(gradv(u)+trans(gradv(u)) );
        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmav = (-idv(p)*Id + 2*M_mu*defv);

        //element_fluid_type Stress(M_Xh,"Stress");
        element_fluid_ptrtype Stress(new element_fluid_type (M_Xh,"Stress"));
        element_fluid_0_type uStress = Stress->template element<0>();

        auto clDef = MECAFLUIDE_CL(M_time);

        ForEachCL( clDef, cl::paroi_mobile,
                   modifVec222(uStress, Stress, Sigmav*N() , PhysicalName ); );

        //modifVec222(uStress, Stress, vec(cst(0.),cst(0.)) , "Entree" );// );
                   //modifVec222(uStress, Stress, -Sigmav*N() , "ParoiH" );
        //uStress.add(1., vf::project( M_Xh->template functionSpace<0>(), markedfaces(M_Xh->mesh(),"ParoiB" ), Sigmav*N() ) );
        //uStress.add(1., vf::project( M_Xh->template functionSpace<0>(), markedfaces(M_Xh->mesh(),"ParoiH" ), Sigmav*N() ) );

        M_Ufluid_ref->container() = Stress->container();
        //*M_Ufluid_ref = *Stress;
        //M_Ufluid_ref->zero();M_Ufluid_ref->add(1.,*Stress);

#else
        auto solFluid = this->getVelocityInRef();
        element_fluid_0_type u = solFluid->template element<0>();
        element_fluid_1_type p = solFluid->template element<1>();
        //Tenseur des deformations
        auto defv = 0.5*(gradv(u)+trans(gradv(u)) );
        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmav = (-idv(p)*Id + 2*M_mu*defv);

        //verifier M_aux est defini sur le ref
        auto Fa = Id+gradv(*M_aux);
        auto Fa11 = trans(Fa*oneX())*oneX();
        auto Fa12 = trans(Fa*oneY())*oneX();
        auto Fa21 = trans(Fa*oneX())*oneY();
        auto Fa22 = trans(Fa*oneY())*oneY();
        //sans le determinant devant car il s annule avec un autre terme
        auto InvFa = mat<2,2>( Fa22,-Fa12,-Fa21,Fa11);


        element_fluid_ptrtype Stress(new element_fluid_type (M_Xh_ref,"Stress"));
        element_fluid_0_type uStress = Stress->template element<0>();
        auto clDef = MECAFLUIDE_CL(M_time);
        ForEachCL( clDef, cl::paroi_mobile,
                   modifVec222(uStress, Stress, Sigmav*trans(InvFa)*N() , PhysicalName ); );


#endif


        if (M_LookComment) std::cout << "[NavierStokes] : getNormalStress finish\n";

        return Stress;//->template element<0>();
    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    typename NavierStokes<Dim, Order, Entity>::element_fluid_ptrtype
    NavierStokes<Dim, Order, Entity>::getVelocityInRef()
    {
        auto solFluid = getSolution();

        //element_fluid_0_type u = solFluid->template element<0>();

        element_fluid_ptrtype Velo(new element_fluid_type (M_Xh_ref,"Velo"));

        Velo->container()=solFluid->container();

        return Velo;


    }



    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::runPtFixe()
    {

        using namespace Feel::vf;

        std::cout << "ptFixe start\n";

        //----------------------------------------------------------------------------------------------------------//

        if ( this->application()->vm().count( "help" ) )
            {
                std::cout << this->application()->optionsDescription() << "\n";
                return;
            }

        //----------------------------------------------------------------------------------------------------------//

#if 0
        Teuchos::ParameterList My_List;
        My_List.set("Max Iters", 1550);
        My_List.set("Tolerance", 1e-10);
        My_List.set("Solver", "GMRES");

        Teuchos::ParameterList& Prec_List = My_List.sublist("Preconditioner");
        Prec_List.set("Type", "ILU");
        Prec_List.set("Drop Tolerance", 1e-3);
        cout<< "\nOLA" <<My_List << endl;

        //( static_cast<BackendTrilinos>(*M_backend)).set_options(My_List);
        BackendTrilinos bbb;
        bbb.set_options(My_List);
#endif

        double erreur=1;
        double err_tol=1e-6;//this->vm()["tolerP"].template as<double>();
        double err1,err2,norm1,norm2;
        uint NbMaxItPtFixe=30;
        uint cptItPtFixe=0;
        uint cptTime=0;


        //U_fluid->container() = M_bdf_fluid->unknown(0).container();

        element_fluid_ptrtype Uoldb( new element_fluid_type(M_Xh,"uoldb"));
        erreur=1;
        cptItPtFixe=0;

        vector_ptrtype R( M_backend->newVector( M_Xh ) );
        vector_ptrtype T( M_backend->newVector( M_Xh ) );
        vector_ptrtype S( M_backend->newVector( M_Xh ) );


        vector_ptrtype F( M_backend->newVector( M_Xh ) );
        sparse_matrix_ptrtype A( M_backend->newMatrix(M_Xh,M_Xh) );

        element_fluid_velocity_type uYYY(M_Xh->template functionSpace<0>(),"xxx");

        while (erreur>=err_tol && cptItPtFixe < NbMaxItPtFixe)
            {
                *Uoldb = *U_fluid;
                //std::cout<<"\nLe Uold Norm L2 : " << Uoldb->l2Norm() << "\n";
                //this->updatePtFixe(*U_fluid,*Uoldb,*M_mesh_velocity,*M_aux);

                this->updatePtFixe(*Uoldb,A,F);
                //this->updateStokes(*Uoldb,A,F);

                std::cout << "PtFixe : Solve start\n";
                vector_ptrtype X( M_backend->newVector( M_Xh ) );
                M_backend->solve( _matrix=A, _solution=X, _rhs=F/*,_maxit=23456*/ );
                //bbb.solve( A,A,X,F/*,_maxit=23456*/ );
                *U_fluid = *X;
                std::cout << "PtFixe : Solve finish\n";


                *T = *U_fluid;
                //this->updateResidual( T, R );
                //modifVec222(uYYY,R, vec(cst(0.),cst(0.)),"Entree");
                //modifVec222(uYYY,R, vec(cst(0.),cst(0.)),"ParoiH");
                //modifVec222(uYYY,R, vec(cst(0.),cst(0.)),"ParoiB");
                //std::cout<<"\nLe residu du ptfixe Norm L2 : " << R->l2Norm() << "\n";

                erreur = /*R->l2Norm();*/this->diffL2(*U_fluid,*Uoldb);
                std::cout << "Erreur : " << erreur <<"\n";
                ++cptItPtFixe;

            }
        std::cout << "ptFixe finish\n";
    }



    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::addLagrangeMult(mpl::bool_<true> /**/, sparse_matrix_ptrtype & A )
    {

        using namespace Feel::vf;

        this->changeRepository();

        auto mesh = this->mesh();
        auto Xh = this->functionSpace();

        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = V.template element<0>();
        element_fluid_1_type p = V.template element<1>();
        //element_fluid_2_type lambda = getLagMult(mpl::bool_<useLagMult>(), V);//V.template element<2>();
        element_fluid_2_type lambda = V.template element< mpl::if_<mpl::bool_<useLagMult>,mpl::int_<2>,mpl::int_<1> >::type::value >();

        form2( Xh, Xh, A ) +=
            integrate( elements(mesh), id(p)*idt(lambda) + idt(p)*id(lambda) );

    }

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::addLagrangeMult(mpl::bool_<false> /**/, sparse_matrix_ptrtype & A )
    {

    }


    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateCLDirichlet(vector_ptrtype& U)
    {

        auto condlim = MECAFLUIDE_CL(M_time);
        auto dirElem = this->getDirichletElement();

        ForEachCL( condlim,cl::dirichlet_vec,
                   modifVec222(*dirElem,U, Expression , PhysicalName ) );
        ForEachCL( condlim,cl::paroi_mobile,
                   modifVec222(*dirElem,U, vf::idv(this->meshVelocity2()) , PhysicalName ) );

        //ForEachCL( condlim,cl::paroi_mobile,std::cout<<"\nololol"<<PhysicalName<<"\n"; );
    }


} // end namespace NavierStokes

} // end namespace Feel

// instantiation
template class Feel::FSI::NavierStokes<MECAFLUIDE_DIM,MECAFLUIDE_ORDER_VELOCITY,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;

//template class Feel::NavierStokes<3,2,Feel::Simplex>;

